Дан график потребительских расходов в Австралии за период от 1989 до 1995 года по кварталам в единицах миллионов долларов. Требуется оценить дальнейшее поведение графика.

Загружаем данные и визуалилируем их:

Попробуем поделить данные на число дней в квартале с целью устранения шума

Ничего не изменилось, поэтому будем рассматривать исходные данные, без деления.

Построим STL-декомпозицию

Значения ошибок такие же больше по амплитуде, как и сезонность. Это вызвано тем, что дисперсия графика растет со временем, что хорошо видно на исходных данных. Чтобы исключить этот эффект, нужно применить преобразование Бокса-Кокса (и визуалилировать результат преобразования).

## [1] 0.1289189

Теперь дисперсия не увеличивается с изменением времени.

ARIMA

Ручной подбор модели

Сделаем ручной подбор модели ARIMA, для этого нужно сначала проверить стационарность данных

print(kpss.test(box_cox_transformed_series)$p.value)
## Warning in kpss.test(box_cox_transformed_series): p-value smaller than
## printed p-value
## [1] 0.01

Получаем, что ряд нестационарен. Применим сезонное дифференцирование:

Проверим снова на стационарность:

print(kpss.test(diff_series)$p.value)
## [1] 0.06700494

Конечно, p получилось больше 0.05, но все равно не сильно. На всякий случай лучше еще раз продифференцировать.

И снова проверяем стационарность.

## Warning in kpss.test(double_diff_series): p-value greater than printed p-
## value
## [1] 0.1

На этот раз p большое, можно считать, что ряд стационарен.

Построим ACF и PACF

Исходя из графиков, можно рассмотреть в качестве начального приближения модели ARIMA(1,1,1)(2,1,1)\(_4\). Рассмотрим несколько близких к ней моделей и посчитаем AICc для каждой:

Модель AICc
ARIMA(1,1,1)(2,1,1)\(_{4}\) -491.3076563
ARIMA(0,1,1)(2,1,1)\(_{4}\) -490.0088317
ARIMA(1,1,0)(2,1,1)\(_{4}\) -489.9499495
ARIMA(1,1,1)(2,1,0)\(_{4}\) -479.9596413
ARIMA(1,1,1)(1,1,1)\(_{4}\) -491.3678759
ARIMA(1,1,2)(2,1,1)\(_{4}\) -489.13115
ARIMA(1,1,1)(2,1,2)\(_{4}\) -490.3129881
ARIMA(1,1,1)(3,1,1)\(_{4}\) -492.127228
ARIMA(0,1,0)(2,1,1)\(_{4}\) -491.4155876
ARIMA(1,1,1)(1,1,0)\(_{4}\) -459.3412136
ARIMA(1,1,0)(1,1,1)\(_{4}\) -490.0123476
ARIMA(0,1,1)(1,1,1)\(_{4}\) -490.0526418
ARIMA(1,1,1)(0,1,1)\(_{4}\) -493.3923725
ARIMA(2,1,1)(0,1,1)\(_{4}\) -491.2466169
ARIMA(1,1,2)(0,1,1)\(_{4}\) -491.2456675
ARIMA(1,1,1)(0,1,2)\(_{4}\) -491.4185055

Самая лучшая модель – ARIMA(1,1,1)(0,1,1)\(_4\)

Построим графики остатков

Построим достигаемые уровни значимости критерия Льюнга-Бокса для остатков

Тут все хорошо, теперь построим QQ-график и гистограмму остатков

Остатки похожи на нормальные. Но это, вместе со стационарностью и несмещенностью можно проверить с помощью критериев:

Гипотеза Критерий Результат проверки Достигаемый уровень значимости
Нормальность Шапиро-Уилка не отвергается 0.0590203
Несмещённость Уилкоксона не отвергается 0.849094
Стационарность KPSS не отвергается 0.1

Разобъем нашу выборку на обучающую и тестовую и посмотрим результат работы прогнозатора:

##                       ME      RMSE       MAE         MPE      MAPE
## Training set   -17.17858  360.8368  263.0813 -0.04074299 0.7666702
## Test set     -1006.41996 1144.9687 1006.4200 -1.76829683 1.7682968
##                   MASE       ACF1 Theil's U
## Training set 0.2143337 -0.1205521        NA
## Test set     0.8199355  0.4816714 0.3549949

Хороший результат.

Автоматическая модель ARIMA

Обучаем автоматическую ARIMA и строим графики остатков

## Series: time_series 
## ARIMA(1,0,0)(2,1,1)[4] with drift         
## Box Cox transformation: lambda= 0.1289189 
## 
## Coefficients:
##          ar1     sar1     sar2     sma1   drift
##       0.9127  -0.0594  -0.1584  -0.6406  0.0353
## s.e.  0.0402   0.1251   0.1084   0.1128  0.0027
## 
## sigma^2 estimated as 0.001566:  log likelihood=254.78
## AIC=-497.57   AICc=-496.94   BIC=-479.92

AICc слегка лучше, чем в ручной модели.

Построим достигаемые уровни значимости критерия Льюнга-Бокса для остатков

Тут все хорошо. Построим QQ-график и гистограмму остатков

Остатки похожи на нормальные. Проверим нормальность, несмещенность и стационарность критериями:

Гипотеза Критерий Результат проверки Достигаемый уровень значимости
Нормальность Шапиро-Уилка не отвергается 0.1601078
Несмещённость Уилкоксона не отвергается 0.8359361
Стационарность KPSS не отвергается 0.1

Все хорошо.

Разобъем выборку на обучающую и тестовую и посмотрим качество прогноза

##                       ME      RMSE       MAE         MPE      MAPE
## Training set   -17.17858  360.8368  263.0813 -0.04074299 0.7666702
## Test set     -1006.41996 1144.9687 1006.4200 -1.76829683 1.7682968
##                   MASE       ACF1 Theil's U
## Training set 0.2143337 -0.1205521        NA
## Test set     0.8199355  0.4816714 0.3549949

Прогноз хороший.

Сравнение ручной аримы и автоматической

Для ставнения двух полученных моделей сначала построим взаимный график их остатков

Теперь прогоним критерий Диболда-Мариано

## 
##  Diebold-Mariano Test
## 
## data:  resres.auto
## DM = 0.43677, Forecast horizon = 1, Loss function power = 2,
## p-value = 0.663
## alternative hypothesis: two.sided

Критерий не говорит нам о том, что одна модель существенно лучше другой. Но, посмотря на AICc, на нормальность гистограмм остатков, видим, что все-таки автоматическая модель лучше.

Прогноз ETS

Построим прогноз ETS

ets_model = ets(time_series, lambda=box_cox_lambda)
print(ets_model)
## ETS(A,A,A) 
## 
## Call:
##  ets(y = time_series, lambda = box_cox_lambda) 
## 
##   Box-Cox transformation: lambda= 0.1289 
## 
##   Smoothing parameters:
##     alpha = 0.7059 
##     beta  = 1e-04 
##     gamma = 0.1704 
## 
##   Initial states:
##     l = 19.4323 
##     b = 0.0351 
##     s=-0.0182 -0.1311 0.1987 -0.0494
## 
##   sigma:  0.0394
## 
##       AIC      AICc       BIC 
## -197.4138 -196.0705 -170.6855

Посмотрим на остатки:

Построим достигаемые уровни значимости критерия Льюнга-Бокса для остатков

Картина выгядит хуже, чем для аримы, но приемлима.

Построим QQ-график и гистограмму остатков

Остатки выглядят нормальными. Проверим нормальность, несмещенность и стационарность критериями:

Гипотеза Критерий Результат проверки Достигаемый уровень значимости
Нормальность Шапиро-Уилка не отвергается 0.119335
Несмещённость Уилкоксона не отвергается 0.808132
Стационарность KPSS не отвергается 0.1

Все хорошо.

Попробуем разбить временной ряд на обучающий и тестовый и посмотреть качество прогноза

##                       ME      RMSE       MAE         MPE      MAPE
## Training set   -18.69817  362.1703  273.7808 -0.01541468 0.8115788
## Test set     -1660.24570 1733.9249 1660.2457 -2.86729008 2.8672901
##                   MASE       ACF1 Theil's U
## Training set 0.2230506 0.04478185        NA
## Test set     1.3526106 0.45730178 0.5329765

Хороший прогноз.

Сравнение ETS и Arima

Сначала построим взаимный график остатков

Теперь попытаемся прогнать критерий Диболда-Мариано

## 
##  Diebold-Mariano Test
## 
## data:  res.autores.ets
## DM = -0.27691, Forecast horizon = 1, Loss function power = 2,
## p-value = 0.7823
## alternative hypothesis: two.sided

Он не обнаруживает качественных отличий, однако, достигаемые уровни значимости критерия Льюнга-Бокса для остатков были лучше у аримы, и показатель AICc тоже. И, наконец, при разбиении на обучающую и тестовую временную серию, арима дала лучший результат.

Поэтому выбираем в качестве лучшей модели автоматическую ариму.

Прогнозирование

##         Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 1995 Q4       63897.38 63104.15 64613.07 62455.42 65094.46
## 1996 Q1       67588.68 66495.47 68630.62 65899.40 69127.77
## 1996 Q2       62118.28 60966.45 63241.71 60216.40 63788.01
## 1996 Q3       64245.70 62952.68 65548.64 62171.56 66158.15
## 1996 Q4       65856.31 64275.27 67390.14 63380.12 68156.67
## 1997 Q1       69765.51 67939.59 71541.78 66894.58 72461.86
## 1997 Q2       64141.11 62346.52 65920.73 61301.19 66783.56
## 1997 Q3       66289.89 64341.89 68209.24 63323.91 69189.15
## 1997 Q4       68108.26 66013.48 70212.74 64816.16 71234.89
## 1998 Q1       72137.36 69799.97 74451.22 68504.67 75516.85
## 1998 Q2       66298.53 64055.68 68549.32 62917.96 69552.63
## 1998 Q3       68575.01 66174.31 70850.28 64954.50 72037.66
## 1998 Q4       70472.80 67902.15 72953.37 66558.60 74192.70
## 1999 Q1       74602.43 71829.34 77262.34 70258.35 78693.10
## 1999 Q2       68591.96 65900.90 71130.89 64555.09 72579.67
## 1999 Q3       70939.11 68075.44 73607.25 66666.48 75207.75
## 1999 Q4       72864.54 69871.68 75737.11 68320.76 77362.23
## 2000 Q1       77115.95 73868.06 80213.78 72110.21 81889.19
## 2000 Q2       70936.17 67849.12 73917.26 66180.91 75477.04
## 2000 Q3       73341.71 70082.74 76516.40 68406.04 78207.06
## 2000 Q4       75321.19 71856.37 78668.49 70092.75 80582.64
## 2001 Q1       79700.07 76022.15 83288.75 74064.07 85142.69
## 2001 Q2       73338.30 69874.32 76725.75 67981.31 78565.90
## 2001 Q3       75814.70 72138.40 79340.52 70319.03 81237.76
## 2001 Q4       77857.07 74003.56 81628.14 72066.86 83703.78
## 2002 Q1       82363.27 78253.63 86398.46 76246.74 88526.70
## 2002 Q2       75814.44 71896.76 79632.21 69944.86 81576.43
## 2002 Q3       78365.17 74316.33 82315.97 72228.45 84464.80
## 2002 Q4       80467.22 76231.64 84576.27 74146.77 86890.15
## 2003 Q1       85103.66 80615.51 89450.63 78331.31 91909.59